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1.  INTRODUCTION 

The  purpose  of  this  study  Is  to  Investigate  a  use  of  time 
series  In  Improving  weather  forecasting. 

Let  us  consider  a  simple  model.  Let  us  assume  that  the 
relevant  variables,  or  functions,  such  as,  say,  the  500  mb  height 
field,  can  be  expressed  In  terms  of  a  set  of  orthogonal  basis 
functions,  whose  coefficients  define  the  field.  There  Is  also  an 
associated  numerical  Integration  program  which  yields  predictions 
of  this  variable  In  the  form  of  predicted  coefficients.  This  method 
of  forecasting  Is  assumed  to  be  more  reliable  than  any  regression 
scheme.  At  regular  Intervals  some  observations  of  the  variable 
are  made  and  the  predicted  coefficients  are  then  updated,  or 
corrected . 

These  corrections  to  the  coefficients  In  the  forecast  are 
analyzed  for  trends.  The  corrections  for  each  (complex)  coeffi¬ 
cient  of  a  (complex)  eigenfunction  or  basis  function  are  analyzed 
separately.  The  purpose  Is  to  predict  the  corrections. 

The  basic  approach  Is  to  define  a  best-fitting  difference 
equation  and  a  best  set  of  Initial  conditions,  by  minimizing 
weighted  sums  of  squares  of  the  residuals,  and  then  use  the 
difference  equation  and  the  Initial  values  to  predict  the  correc¬ 
tion.  These  are  used  In  a  secondary  routine  to  Improve  the  forecast. 

Testing  of  the  method  In  the  500  mb  field  with  a  simple 
non-dl vergent  model  due  to  Bourke  (1972)  shows  skill  with  some 
coefficients . 
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OUTLINE  OF  THE  PROBLEM 

Let  us  assume  that  we  have  some  function,  or  functions,  which 
e  representing  In  terns  of  an  orthonormal  set  of  basis 
functions.  In  particular,  let  us  consider  the  500  mb  height 
function,  or  perhaps,  the  stream  function,  expressed  on  the  surface 
of  the  earth  In  terms  of  spherical  harmonics.  These  are  expressed 
as  complex  trigonometric  functions  of  longitude  and  Legendre 
functions  of  the  sine  of  the  latitude.  At  time  tN  the  coefficients 
of  the  spherical  harmonics  form  a  vector  or  column  matrix  ;  Its 
components  are  complex  numbers  ak  N,  k«l ,  ....  K,  If  the  series 
Is  truncated  to  K  terms.  The  vectors  AN,  and  each  of  the  components 
ak,N*  forM>  a  *1me  ser1es»  or  sequence. 

The  sequence  Is  generated  conceptually  as  follows.  Let  us 
assume  that  at  time  t^j  we  have  the  matrix  An_j  of  coefficients 
for  the  analyzed  field.  This  Is  considered  to  represent  the  true 
field  at  that  time.  We  have  also  an  Integration  procedure  by  which 
we  can  generate  predictions  A^  N+1 .  ....  for  the  future 

values  Aj|,  AN+1 ,  ....  as  many  as  we  feel  are  useful  or  reliable. 

Then  at  time  tN,  12  hours  later,  observations  are  made.  These  are 
combined  with  the  predictions  A^_^  N  toylelda  new  analyzed  field, 
and  a  new  set  of  coefficients  A^.  The  corrections  (changes, 
errors,  or  discrepancies). 

Ip  *  An  ~  Ap_^  n  •  2,  3,  ...,  N,  (2.1) 

also  define  a  time  series  or  a  sequence.  If  we  could  In  some  way 
analyze  this  sequence,  Z2,  Z3,  ..,  ZN,  and  from  It  predict 
Zw+1 •  ZN+2,  ....  ZN+m  reliably  then  we  could  predict  ahead  m  steps 


_ 
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reliably.  Our  purpose  Is  to  predict  mo  re  accurately,  or  More 
reliably,  and  perhaps  to  extend  the  range  of  prediction.  We  shall 
see  that  there  are  two  ways,  at  least,  that  this  Improvement  might 
be  effected.  Both  depend  on  the  use  of  time  series  as  developed 
In  the  next  sections  . 

The  Z's  and/or  their  components  form  the  data  for  the  analysis. 
As  such  they  are  called  observations  (not  to  be  confused  with 
meteorological  observations).  The  terms  component  and  coefficient 
are  used  Interchangeably.  Vectors  have  components  and  basis 
functions  have  coefficients  when  they  are  used  In  the  representation 
of  another  function.  The  components  afc  n  of  A^,  for  example,  are 
the  coefficients  In  the  representation  of  some  meteorological 
variable.  Basis  functions  are  often  called  modes,  particularly  in 
mechanics . 


3.  INTRODUCTORY  TINE  SERIES,  SIMPLIFYING  ASSUMPTIONS 


Let  us  consider  the  sequence  of  corrections,  Zn>  n  *  2,  3,  ..., 
The  vector  ZR  has  K  complex  components.  There  is  one  complex  number 
zk  n  ^or  each  the  K  complex-valued  spherical  harmonics  that  Is 
used.  To  save  writing,  let  z  .  rather  than  z.  .  denote  a  typical 

n  K  1 0 

coefficient;  we  will  analyze  each  of  these  sequences  individually. 
Part  of  the  rationale  is  that  in  a  linear  uncoupled  system  the 
various  modes  move  Independently;  In  a  loosely  coupled  system  we 

I  *  *  •  i  I  !  I 

would  expect  similar  behavior.  There  are  other  Important  reasons 

which  we  will  take  up  later.  We  shall  usually  consider  z_  to  be  a 

n 

pair  of  real  numbers,  a  two-by-one  matrix. 


Let  us  now  assume  that  zn  can  be  expressed  as  a  sum  of  two 
terms,  described  below. 


Here  yn  Is  a  random  variable.  Its  components  each  are  assumed  to 

2 

have  expected  value  zero  and  expected  square  o  ,  Independent  of  n; 
they  are  also  assumed  to  be  uncorrelated.  The  vector  xn  Is  assumed 
to  be  a  solution  to  a  difference  equation  of  the  form 


The  coefficients  Cj ,  ....  Cp  are  unknown  two-by-two  matrices,  and 
the  order  p  Is  to  be  chosen  some  way.  For  the  present  let  us  take 
p  ■  2,  which  seems  at  this  time  to  be  a  likely  value.  The  C's  are 
not  constants  necessarily,  but  may  vary  adaptively  with  N. 

To  determine  the  coefficients  In  (3.2),  let  us  choose  the  C's 
to  minimize  a  residual  of  the  form 


where  the  wn  are  positive  weight  factors  to  be  chosen  or  determined; 
we  will  discuss  this  later.  To  minimize  R  let  us  differentiate  with 


respect  to  the  elements  of  the  C's;  we  get  a  system  of  equations  of 


the  prime  (')  Indicates  the  transpose  of  a  matrix.  The  terms  like 
*_  «  represent  two-by-two  matrices  so  that  the  elements  of  the 

n  -  i  n  -  p 

Matrices  are  themselves  matrices  of  order  two. 

Another  problem  Is  the  choice  of  the  weight  factors  wR.  For 

a  time-independent  process  we  should  use  equal  weightings,  which 

simplifies  several  relations.  However  the  fundamental  relations 

may  vary  rather  rapidly  with  time  sometimes,  and  these  are  the  very 

times  that  are  most  critical.  One  way  to  take  account  of  this  Is 

to  use  exponential  weightings  In  (3.3)  and  (3.4).  After  some 

starting  routine  we  weight  each  new  observation  with  a  uniform 

value  w  and  decrease  all  earlier  weightings  then  by  a  factor  1-w. 

Me  do  this  as  follows.  Let  QN  be  the  matrix  of  coefficients  In 

(3.4);  It  has  two-by-two  block  elements  of  the  form 
N 

’-J'1 . pi  (3-5> 

the  right-hand  side  of  (3.4)  has  elements  Q  When  we  get  the 
observation  z^,  we  update  j  N  -|  by  the  relation 

Q1,j,M  ’  Q1.j,N-l  +  w(zN-1zN-j  '  Q1,j,N-l)-  (3‘6) 

Me  may  actually  update  Q  by  first  shifting  each  block  element 
down  and  to  the  right  one  place,  dropping  out  the  last  row  and  the 
last  column.  Then  we  update  the  first  (block)  row,  using  (3.6), 
and  then  the  first  column  by  symmetry.  The  weight  of  an  observa¬ 
tion  from  time  t  for  the  coefficients  calculated  at  time  tM  Is 

n  n 

w(l-w)N‘n;  we  see  that  weightings  decrease  exponentially  with  time. 
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The  best  choice  of  w  Is  another  problem.  If  w  Is  too  small 
the  system  Is  slow  to  respond  to  changes,  and  If  w  Is  too  large, 
random  errors  cause  excessive  errors  In  the  predicted  values. 

The  above  procedure  for  updating  the  coefficients  In  (3.2)  1$ 
rather  efficient,  both  In  terms  of  storage  and  computation.  The 
matrix  Q  Is  always  singular  until  at  least  2p  observations  have 
been  made.  It  Is  usually  positive  definite  after  that,  but  It  was 
sometimes  found  to  be  near  singular,  especially  for  the  smaller 
values  of  N. 

In  order  to  predict  the  z's  ahead,  using  (3.2),  we  also  need  p 
starting  or  Initial  values;  we  can  take  these  conceptional ly  as 

estimates  of  xN,  xn_j . xN-p+l  *  The  determ1 nati on  °f  these  Is 

taken  up  In  the  next  section;  It  requires  more  storage  and  more 
computation  than  the  above. 

In  the  following  section  we  will  generally  consider  the  case 
p  *  2. 

There  may  be  some  minor  discrepancies  in  the  first  Index;  for 
example,  when  the  a's  are  Involved  we  tend  to  think  of  z2  as  the 
first  of  the  z's,  and  later,  when  treating  the  z's  we  think  of  the 
first  as  z^  . 


We  can  simplify  the  programming  somewhat  If  we  define  the 
matrices 


D"  ■  (*„.»„) 

U  *  ( x^ ,  x^  )  1  *  ( xi  2  •  x22  *  ^11*  X21  ^  '  • 


(4.3) 


Is  the  1'th  component  of  Xj.  Then  Dn  satisfies  Eq.  (3.2)  and 


has  the  Initial  values 
.n 


C^Dn”^  ♦ 


D1  -  (0,  I),  D2  »  (I,  0) 


(4.4) 


Equation  (4.1)  now  reduces  to 


zn  -  D  U 


(4.5) 


If  we  knew  the  four  components  of  U,  we  could  get  starting  values 
*N*  *N-1  from  Eq  *  (3  -2)  • 

xn  =  C1 xn- 1  +  C2xn- 2  *  (4,6) 

and  If  we  select  or  estimate  U  in  some  way  we  get  corresponding 
estimates  for  xN,  x^  j  .  Of  course  It  would  be  nicer  to  determine 
these  in  a  more  direct  way,  but  there  seems  to  be  no  way  to  do 
this  . 

Now,  by  assumption,  the  expected  value,  E  y^  =0.  Hence  it 
seems  reasonable  that  the  best  choice  for  U  Is  one  which  minimizes 
a  weighted  sum  of  squares  for  the  estimates  of  yi ,  say. 


,  N  ,  ,  N  _  2 

R  ‘  2  Z  Vn  =  2  Z  wn(D  U  ’  zn> 


(4.7) 


To  effect  this  we  differentiate  with  respect  to  the  components  of 
U;  If  we  call  Its  components  u j ,  the  j'th  partial  derivative  yields 

Rj  3  z  wn[d1,j^  d1,kuk  '  zl,n)  +  d2,j(^  d2,kuk  '  z2,n^ 

-  0,  o  (4.8) 


where  djk  ere  the  components  Db.  We  solve  these  equations  for  U, 
which  defines  Xj,^,  and  from  these  we  get  our  best  estimates, 
xN_l  ,  iN  for  starting. 

There  are  several  problems.  First,  the  sequence  zR, 
n  3  1,  2  ...,  N,  may  be  very  long,  so  that  computing  times  and 
errors  may  be  significant.  Further  the  observations  z„  for  small 
values  of  n  may  not  be  relevant  to  today's  weather.  We  also  need 
to  recall  many  observations.  These  problems  suggest  limiting  the 
range  of  summation  In  (4.7)  to,  say,  the  last  ten  terms,  which  we 
did. 

Another  problem  needs  some  explanation.  We  really  just  need 

xN1  ,  xN  to  predict  ahead.  However  we  cannot  get  these  directly. 

It  really  does  not  matter  whether  we  find  xN and  xN,  or  xN_g 

and  x^g;  we  get  the  former  from  the  latter  by  use  of  difference 

Eq.  (3.2).  We  might  perhaps  try  to  rewrite  the  difference 

equation,  to  solve  for  x„  ,  in  terms  of  x„  ,  and  x  .  to  avoid 

n -  c  n- 1  n 

this.  However  we  must  invert  the  matrix  C2  to  effect  this,  and 
there  Is  no  reason  that  it  cannot  be  singular,  or  nearly  so, 
randomly . 

In  the  program  we  have  15  %’s  actually  Indexed  from  1  to  15. 
The  first  10  of  these  are  the  10  approxlmants  to  the  last  10 
observations  zN_g,  •••»  z^;  the  last  five  are  predictions  for  the 
future  values  z^+1 ,  ...»  zN+5. 
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5.  TWO  SU66ESTED  METHOOS  FOR  IMPROVING  PREDICTION  BY  THE  USE  OF 
A  TINE  SERIES 

There  are  two  particular  methods  which  seem  feasible  for 
Improving  the  forecast.  Each  of  these  seems  feasible,  but  each 
poses  some  problems. 

5.1  Method  1 

Let  us  consider  again  the  matrices  for  the  coefficients  In  the 

analyzed  fields,  An»  n»l .  N;  we  wish  to  Improve  the  predictions 

of  these.  Let  us  assume  that  we  have  generated  and  stored  the 
following.  First,  we  have  stored  the  last  ten  corrections 
*N-9’  **•*  *n*  ^ave  a^so  stored,  for  each  term  In  the  Z's, 
the  matrix  Q  which  defines  the  coefficients  Cj  In  the  difference 
equation . 

The  routine  Is  the  following  for  each  component  of  Z.  As 
soon  as  the  new  value  of  Z,  ZN,  Is  obtained,  we  store  It  and 

discard  the  old  Z^g.  Then  we  update  the  matrix  Q  and  solve  for 

the  coefficients  of  the  difference  equation  as  described  In 
Section  3.  Then  we  solve  for  the  Initial  values  xN-1 »  xN,  as 
described  In  Sect.  4,  and  finally  we  predict  ahead  to  get 

XN+1 *  xN+2*  **■*  *N+M*  as  raany  as  **  wish.  These  define  the 
x's ,  or  X's . 

Next  let  us  use  the  dynamic  equation  and  Integrate  ahead  one 
time  step,  to  get  an  estimate  AN >f)+1  ,  of  AN+1 .  Then  let,  say. 

*N,N+1  *  AN,N+1  » 
and 

*  ***  A 

AN,N+1  "  AN,N+1  +  XN+1  •  l5*1) 


10 


A 

Me  get  successive  values  In  a  similar  way.  Mhen  we  heve  Ag  g+j 
we  Integrate  to  obtain  Ag tu+j+^ »  •  first  estimate  of  Ag  g+ ^ . 

A 

Then  we  add  Xg+j+1 ,  to  get 

A  w  A 

AN,M+j+l  *  AN,M+j+l  +  XM+j+l •  t5-2) 

A 

Mlth  the  Initial  value  Ag  g  *  Ag,  the  sequence  of  predictions, 
or  forecasts.  Is  defined. 

I  i 

In  the  Initial  part  of  the  study,  a  simple  Integration 
routine  was  used  for  predicting,  Ag  N+g,  a  model  using  the  non- 
dlvergent  barotroplc  vortlclty  equation.  Mhen  this  model  Is  used, 
the  method  of  correcting  above,  alternately  adjusting  and 
Integrating  In  feasible  and  practical. 

Mhen  a  more  accurate  and  complicated  weather  prediction 
model  Is  used,  serious  difficulties  arise  with  this  method. 
Restarting  a  time  Integration  after  rather  arbitrary  adjustments 
to  one  or  more  variables  may  generate  spurious  gravity  waves  which 
degrade  the  prediction  and  offset  the  adjustments.  To  eliminate 
this  difficulty  a  second  method  may  be  used.  In  which  a  separate 
series  Is  generated  for  each  Interval  of  prediction.  The  appro¬ 
priate  term  Is  used  for  each  Interval  and  the  adjusted  functions 
are  never  Integrated. 

5.2  Method  2 

There  Is  another  method  which  gets  around  the  above  problem, 
at  the  expense  of  Increasing  the  storage  and  computation  by  a 
factor  of  roughly  M,  If  we  use  M  different  Intervals  of  prediction. 
That  Is,  If  we  wish,  say  to  adjust  the  12-,  24-,  and  36-hour 
forecasts  we  must,  generate  three  time  series. 

11 


Now  assume  that  we  wish  to  predict  ahead  M  steps.  At  the 
time  tN  we  will  have  the  predictions  based  on  Integration, 


N.N+l  *  *•••  "N.N+M  1  N,N+nrm»l  ,M 
For  the  corrections  we  will  save  M  sequences  Z 


Let  us  consider  any  particular  value  of  m.  To  generate  the  data 
for  the  time  series  we  must  store  each  time  the  uncorrected 


we  will  obtain  and  store 


We  will  thus  have  M  time  series  to  be  analyzed.  Each  Is 
analyzed  as  discussed  earlier.  In  this  case  however  we  will  use 
a  single  prediction  from  each  one:  from  Z™  we  generate  the  single 
correction  for  m  time  steps  ahead.  If  we  denote  the  corrections 
by  X"  the  modified  values  will  be 


This  method  has  the  advantage  that  the  modified  terms  are  not 
Integrated.  It  has  the  disadvantage  that  It  requires  one  time 
series  for  each  value  of  m  used,  which  leads  to  larger  storage 
and  computational  requirements. 


6.  COMMENTS 

The  problem  here  differs  fro*  the  aost  coaaon  applications  of 

tlae  series.  We  have  a  rather  saall  aaount  of  Inforaetlon  which 

/ 

tends  to  be  masked  by  a  large  randoa  eleaent.  Froa  this  we  are 
trying  to  predict  transients,  or  trends,  for  relatively  short 
periods,  perhaps  one  to  ten  tlae  steps.  We  are  not  directly 
Interested  In  the  smoothed  values,  since  we  are  not  concerned  with 
what  has  happened,  except  as  our  ability  to  fit  It  reflects  on  our 
ability  to  predict.  Our  feeling  Is  that  the  variables  we  are 
predicting  may  change  rather  rapidly  so  that  the  weather  two  weeks 
ago  Is  of  little  interest;  even  the  equation  which  governs  the 
behavior  may  well  have  changed.  We  cannot  really  make  use  of 
long-term  observations,  as  for  a  stable  system,  which  allow  more 
accurate  predictions. 

We  are  trying  to  predict  the  short-range  behavior  of  a  part 
of  a  non-linear  system  of  high  order  by  solutions  to  a  linear 
system  of  low  order.  Whether  we  try  to  do  this  by  a  difference 
equation  or  by  Taylor  series  depends  on  the  type  of  data,  and  the 
nature  of  the  expected  solutions.  In  our  case  both  the  form  of 
the  data  and  the  anticipated  periodic  properties  of  the  expected 
solutions  suggest  difference  equations. 

The  papers  by  Jones  (1963),  (1964)  particularly,  and  (1965) 
suggested  the  general  problem  here.  These  treat  more  the 
descriptive  and  the  long-range  forecasting  problem.  The  series 
are  considered  to  be  stationary  In  the  sense  that  all  observations 
are  weighted  equally;  several  formulas  are  then  simplified. 
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The  references  found  most  useful  In  this  study  were  Whittle 
(1963)  end  Gelb  (1974).  Recent  developments  In  the  use  of  spherical 
harmonics  In  meteorology  (see  GARP,  1974)  and  Fast  Fourier  Transforms 
make  the  methods  appear  feasible. 

The  method  Is  clearly  that  of  least  squares.  It  Is  also  a 
linear  regression  method;  all  of  the  equations  to  be  solved  are 
linear.  If  we  were  to  use,  say,  our  estimates  xN  -j  ,  xN  In  a 
routine  to  try  to  Improve  the  C's,  then  It  would  be  non  linear. 

The  solutions  to  the  difference  equations  are  basically  complex 
exponentials,  that  Is,  a  combination  of  real  exponentials  and 
trigonometric  functions. 

There  are  a  number  of  points  that  should  be  discussed  or 
cl  arl  fled . 

6.1  Consistency 

It  should  be  pointed  out  that  even  for  a  stationary  system 
the  solutions  are  not  consistent,  as  follows.  Let  us  assume  that 
we  have  a  solution  to  a  difference  Eq.  (3.1)  to  which  random 
uncorrelated  terms  yn>  each  with  expected  value  0  and  expected 
square  a  Is  added.  Then  Eq.  (3.4)  will  not  yield  the  desired 
coefficients  In  (3.2).  The  expected  elements  of  the  matrix  of 
coefficients  Q  In  (3.5),  from  (3.4)  will  be 
N  2 

"h  “u”  V  '-i-1 . .  (<M) 

here  I2  Is  the  second-order  Identity  matrix,  and  6  Is  the  Kronecker 

6,  (*1  when  1«j  and  0  when  Ipj).  That  Is,  the  Q  matrices  have  an 
2 

extra  term  o  on  the  main  diagonal.  This  suggests  that  If  we  have 
2 

a  measure  of  o  we  could  subtract  It  from  the  elements  on  the  main 
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diagonal  of  Q,  In  a  sort  of  negative  ridge-regression  scheme. 

This  Is  very  risky  procedure  however,  since  It  drives  the  aatrlx 
of  coefficients  toward  singularity.  Difficulty  has  been 
encountered  several  tlaes  because  this  aatrlx  Q  was  singular,  and 
the  danger  Is  aggravated  when  exponential  weighting  with  a  large 
decay  rate  Is  used. 


6 .2  Best  Order  for  the  Difference  Equation 


He  have  suggested  that  we  use  a  second-order  difference 
equation,  p  ■  2.  Soae  of  the  reasons  follow.  Let  us  consider  a 
systea  which  Is  tlae  Independent,  and  assuae  that  have  found  the 
correct  order  p  and  the  true  coefficients  C  In  the  difference 

2 

equations.  In  this  case  the  residuals  yield  an  estlaate  for  a  . 
Now  consider  the  residual 


since  the  x's  satisfy  (3.2).  The  expected  value  of  the  residual 


tv*/  -  °  . 

where  the  last  term  In  the  parenthese  denotes  the  sum  of  the 
squares  of  all  of  the  components  of  all  of  the  C‘s. 

Now  If  we  try  to  fit  the  solutions  to  a  second-order  system 
by  a  third-order  difference  equation.  In  the  absence  of  noise,  the 

since  the  solutions  satisfy 
If  we  have  a  little  noise 


aatrlx  of  coefficients  Is  singular 


third-order  difference  equations 


we  nay  have  a  near-singular  Matrix  and  have  large  values  for  the 
c1.j,k*  1n  th1 s  case  *•  9«t  a  large  expected  value  In  (6.1)  for 
the  residual  .  Limited  numerical  checks  make  the  second-order 
equation  appear  satisfactory. 

6.3  Determining  Initial  Values 


There  was  no  particular  reason  for  choosing  ten  values  to  fit 
when  we  determine  the  starting  values  for  prediction.  It  was  felt 
that  six  would  be  the  smallest  number  to  be  considered,  that  fewer 


would  make  the  methods  vulnerable  to  random  errors  In  the  last 


terms.  Many  more  than  ten  could  cause  roundoff  errors 


6.4  Criteria  for  Adjusting  Terms 


We  do  not  expect  to  adjust  or  correct  all  of  the  coefficients. 
For  some  the  errors  may  be  too  small  to  warrant  the  time  and 
trouble.  Other  coefficients  may  contain  such  random  elements  and 
change  so  unpredlctably  that  trends  are  effectively  obscured.  The 
decision  to  adjust  or  correct  will  probably  depend  on  an  analysis 
of  similar  sequences  In  the  past.  There  are  two  basic  ways  to  do 
this,  one  based  on  a  long  sequence  of  past  data  (that  Is,  on 
climatology),  and  the  other  on  recent  data  (corresponding  to  recent 
weather) . 

We  may  analyze  a  long  sequence  for  a  coefficient,  perhaps  over 
a  period  of  years,  and  from  this  make  a  decision,  beforehand 
always  to  correct,  or  never.  In  the  program. 

An  alternate  way  Is  to  see  how  well  the  smoothed  values  fit 
the  last  ten  corrections.  We  generate  these  values  naturally  In 
determining  the  starting  values  for  prediction,  so  this  entails 
little  computation.  We  may  consider,  for  example,  the  ratio 


10 

1.25  I 
1 


wn<2N-10+n 


,  10 
*  )2/Z 
n  1 


"n2N-10+n 


If  this  Is  Much  less  then  one  we  have  fitted  the  past  ten  values 
well,  and  If  It  exceeds  one  we  have  no  skill  In  fitting  the 
observations.  The  factor  1.25  adjusts  for  the  fact  that  we  could 
fit  two  terms  exactly  with  the  Initial  conditions.  We  might 
correct  or  not  accordingly  as  this  number  was  below  or  above  some 
criterion,  say,  .75.  The  philosophy  Is  that  If  we  can  fit  a 
sequence  of  observations  well,  then  we  can  predict  well.  (More 
logical  Is  the  converse;1f  we  cannot  fit  the  observed  values,  then 
we  cannot  predict  well.)  It  Is  probably  not  worthwhile  to  use 
this  criterion  unless  we  usually  correct. 

6.5  Criteria  for  Improvement 

In  the  numerical  study  we  used  two  principal  criteria  for 
Improvement.  Both  of  these  compared  In  some  way  the  magnitude  of 
the  errors  In  prediction  with  adjustment  with  the  errors  generated 
without  the  adjustment.  In  the  first  criterion,  the  error  In  the 
predicted  values(s)  was  compared  with  the  rms  value  of  the  last 
ten  corrections  (without  adjustment).  This  figure  was  calculated 
as  the  predictions  were  under  way.  For  the  second  criterion,  a 
long  sequence  of  adjustments  were  calculated.  Then  the  rms  values 
of  the  errors  In  prediction  after  adjustment  was  compared  with  the 
rms  value  of  the  error  In  prediction  without  adjustment: 


■I',  - 

the  range  of  summation  of  n  Is  over  each  term  for  which  a  correc¬ 
tion  to  z  was  made.  These  measures  are  easily  made  In  the  machine, 
and  they  define  a  simple  norm  for  the  error  relation. 
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It  Is  not  clear  these  norms  are  good  ones  from  a  meteorological 
point  of  view.  Very  likely  the  ultimate  test  will  rest  on  maps  and 
criteria  for  judging  them. 

6 .6  Real  Coefficients 

Most  of  the  coefficients  are  complex.  The  one  associated  with 
the  modes  that  are  Independent  of  longitude  (merldlanal)  are  real 
and  hence  one  dimensional.  The  routines  developed  then  Involve 
scalar  equations,  rather  than  two-dimensional  vectors.  The 
corresponding  spherical  harmonics  are  the  Legendre  polynomials, 
the  simplest  of  the  Legendre  functions. 

7.  RESULTS 

The  procedure  for  analyzing  the  time  series  was  checked  on 

several  sets  of  made-up  data,  and  on  two  coefficients  ag  and  a4 
12 

(of  Pj  and  P£),  using  a  difference  equation  of  second  order,  p  *  2. 
For  data  made  up  of  solutions  to  a  difference  equation,  with 
random  noise  added,  the  solution  could  be  recovered  If  the  noise 
was  not  too  large,  but  the  routine  did  not  perform  very  well  when 
the  energy  In  the  noise  (rms  value)  was  larger  than  In  the  solution. 

For  meteorological  data,  on  a  string  of  60  pairs  of  correc¬ 
tions,  the  predictions  x^  to  x55  led  to  a  definite  Improvement  for 
one  coefficient.  The  best  value  for  the  weighting  w  of  the  new 
values  seemed  rather  small:  the  values  .025,  .05  and  .1  gave 
similar  results  and  were  better  than  larger  values.  The  numbers 
below  were  determined  as  follows.  He  considered  the  coefficient 
of  p] .  A  set  of  10  data  points  were  used  at  the  start  to  get 
coefficients  for  the  best  difference  equation  and  the  Initial 


values.  Five  values  were  predicted,  corresponding  to  the  12-hour 

correction  to  be  nade  In  12  hours,  24  hours,  ...,  60  hours.  The 

error  In  predicting  these  was  then  calculated  and  normalized  by 

dividing  by  the  rms  of  the  10  original  data.  Then  the  time  Index 

was  advanced  by  one  and  the  operation  repeated. 

The  numbers  In  row  one  are  the  rms  values  of  an+1  -  an  n+1 , 

normalized  by  the  rms  values  of  *n_g»  zn.g»  •  ••»  zn-  The  numbers 

In  row  m  denote  the  rms  value  of  a„._  -  a  normalized  by 

n+m  n+m- i ,n+m 

the  rms  value  of  z„  a,  ....  z  .  Each  Is  then  a  correction  to  be 

n*s  n 

applied  after  Integrating  one  time  step.  Each  entry  comes  from 
the  sum  of  46  terms  . 


No .  of 

Time  Steps 

Weight  w  .025 

.05 

.1 

.2 

.4 

1 

.6239 

.6257 

.6131 

.6404 

.81  38 

2 

.6370 

.6401 

.6312 

.6629 

.8640 

3 

.6670 

.6751 

.6834 

.7599 

1 .3095 

4 

.6573 

.6698 

.6960 

.8314 

1 .5523 

5 

.6513 

.6636 

.6880 

.8496 

1  .9109 

There  are  several  Interesting  points  about  the  data.  First, 
the  results  from  the  small  values  of  w,  from  .025  to  .1  are  all 
comparable.  This  Indicates  a  stable  trend:  when  w  =  .025  the 
first  of  the  10  data  points  had  a  weighting  almost  .8  of  the  last; 
they  are  weighted  nearly  equally.  For  w  =  . l o  the  relative 
weighting  Is  about  .39.  A  surprising  feature  Is  that  all  five  of 
the  predictions  are  so  similar.  This  Is  further  Indication  of  a 
stable  trend,  a  large  "random"  element,  and  a  satisfactory  value 
of  p,  and  suggests  longer  predictions  are  feasible. 
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The  corrections  to  the  fourth  coefficient  Mere  not  predicted 
successfully;  the  error  In  the  predicted  values  was  consistently 
a  little  large  than  the  corrections.  However  various  features  of 
the  solutions  are  quite  similar  to  the  other  solutions;  the 
smaller  values  of  w  give  better  predictions  and  all  five  predicted 
values  are  similar  In  size,  with  later  predictions  often  anomalously 
having  smaller  error  than  the  first. 
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SYMBOLS:  DEFINITIONS.  AND  PLACES  WHERE  INTRODUCED 


Section  2 


*n("ak,n> 


Z«*A  -A„  i  _ 
n  n  n- 1  ,  n 


time  N  Is  usually  associated  with  the  last 
or  next-to-last  observation. 

t  i  ri  X  t  s  *2  Q  .  fjf  3  i>  ij  V  ^  » 

the  matrix,  or  vector,  of  coefficients  for 
the  analyzed  field  at  time  tn. 

the  k'th  element  of  A„,  k»l  ,  ....  K.  Most 

n 

of  the  a.  _  are  complex  numbers  and  are 
k  ,  n 

treated  like  two-dimensional  vectors. 

the  number  of  spherical  harmonics  used  to 
represent  a  function. 

a  typical  element  of  Ar,  written  without 
Index  k  to  save  writing. 

the  value  predicted  for  An^  at  time  tn. 

This  represents  the  discrepancy  between  the 

analyzed  field  at  tfl  and  the  value  predicted 

for  It  at  time  t„  Also  called  a  correction, 

n*  i 


Section  3 


a  typical  element  of  Zn» 

the  component  of  zn  that  Is  assumed  to  be 
random . 

the  component  of  zn  that  Is  assumed  to  satisfy 
a  difference  equation  (see  below). 

two  by  two  matrices  In  the  difference  equation 

n  1  n-1  p  n-p 

weighted  sum  of  residual  errors. 

weighting  of  new  terms. 

the  weighting  of  term  associated  with  tn,  at 

time  N>n;w„  decreases  with  N  for  a  fixed 
—  n  u 

value  of  n:w„  ■  w(l-w) 
n 
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Section  3  (continued) 

xn  a  best  approximation  for  xn,  made  at  some 

definite  time  t^.  (It  might  also  have  been 
designated  by  xk  N  n«) 

X„  the  vector,  or  matrix,  of  x  's. 

n  n 


Section  4 


»*  ■  <VB„> 

X1J 


matrices  defined  In  Eq.  (4.2),  used  In 

fl  ndl  ng  x  ‘  s  . 
n 

(see  Eqs.  (4.3),  (4.4)). 

the  I'th  component  of  x^,  1,j*l,2. 

,x11  ,x]2)  ' 


Section  5 

a 

AN ,N+m 
AN ,N+m 


an  optimal  predictor  for  AN+m, 

a  temporary  estimate  of  AN  N+m 
obtaining 


made  at  tN 
,  used  In 


A 


n-m,n 


typical  elements  of  Z™. 


R 


a  general  symbol  for  a  quadratic  residual 
to  be  ml nlml zed . 
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